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ABSTRACT 

We present a new numerical method of special relativistic resistive magnetohydrodynamics with scalar re- 
sistivity that can treat a range of phenomena, from nonrelativistic to relativistic (shock, contact discontinuity, 
and Alfven wave). The present scheme calculates the numerical flux of fluid by using an approximate Riemann 
solver, and electromagnetic field by using the method of characteristics. Since this scheme uses appropriate 
characteristic velocities, it is capable of accurately solving problems that cannot be approximated as ideal mag- 
netohydrodynamics and whose characteristic velocity is much lower than light velocity. The numerical results 
show that our scheme can solve the above problems as well as nearly ideal MHD problems. Our new scheme 
is particularly well suited to systems with initially weak magnetic field, and mixed phenomena of relativistic 
and non-relativistic velocity; for example, MRI in accretion disk, and super Alfvenic turbulence. 
Subject headings: plasma, relativistic resistive MHD, methods: numerical 



1. INTRODUCTION 

The magnetohydrodynamics (MHD) approximation has 
some interesting properties, for example, the flux freezing 
and magnetic pressure; the former can be used for the col- 
limation of the jet, and the latter for the acceleration of the 
plasma. Thus, the magnetic field is considered an essential 
ingredient for many astrophysical phenomena. In particular, 
many observations indicate that most of the high energy phe- 
nomena in astrophysics are related to the strongly magnetized 
relativistic plasma arou nd sor ne compact obiects, for exam- 
ple, AGN |Antonucci"1993'; lUrry & Padovani'liggS), rela- 
tivistic jet iBlandford & Konigl 1979; Mirabel & Rodriguez 
[1999^, pulsar wind (Rees & Gunn 1974; Camus et al. 2009), 
gamma-ray bursts (Woosley 1993; Piran 2004), and so on. 
Since it is extremely difficult to solve the relativistic MHD 
(RMHD) equations analytically, the theoretical investigations 
in fully nonlinea r regimes are mainly based o n the numer- 
ical s imulations dMcKinnev & Gammie I l2004t llnoue et aUI 
I2OIOI) . Most of these studies approximates the plasma as 
the ideal RMHD fluid. One reason for this is that the ideal 
RMHD is an excellent approximation of high energy phe- 
nomena for ordinary parameters. However, when one con- 
siders extreme phenomena, such as the neutron star mergers, 
or the central engines of GRB, the electrical conductivity can 
be small, and highly resistive regions may appear In addi- 
tion, when one considers the magnetic reconnection, the re- 
sistivity plays an essential role in this phenomenon. Magnetic 
reconnection is one of the most important phenomena, since 
it is highly d ynamic, and it chang es mag netic field energy into 
fluid ener gy (Zweibe l & Yamada 120091: IZenitani et al. 120091: 
iLietal. 112007) . Though numerical results of ideal RMHD 
exhibit magnetic reconnection, this originates in the purely 
numerical resistivity, and this is unphysical. For this reason, 
using resistive RMHD is important for the understanding of 
reconnection and related phenomena. 

In order to consider Ohmic dissipation, one only has to 
take into account an additional term -V x ( V x B)/cr in the 
induction equation of non-relativistic MHD. However, simi- 



lar to other non-relativistic dissipation, this induction equa- 
tion is parabolic and it is well-known that this equation is 
acausal. As a result, if one takes into account Ohmic dissi- 
pation in a relativistic MHD in a similar way, the equation 
inevitably includes unphysical exponential growing modes, 
and uns table for small perturbations simi lar to other dissi- 
pation dffiscock & Lindblomlll983L Il985h . This unphysical 
divergence results from the fact that one neglects the time 
derivative of the electric field in the induction equation with 
Ohmic dissipation. For this reason, when one takes into ac- 
count the Ohmic dissipation, one has to consider the time 
evolution of the electric field, that is, one has to deal with 
the relativistic electromagnetic hydrodynamic equation. This 
equation is a telegrapher equation, and satisfies the causality. 

In this paper, we present a new numerical scheme for the 
resistive RMHD. There are several examples of pioneering 
work for resistive RMHD, for example, Komissarov (2007, 
hereafter K07) proposed numerical method that solves hyper- 
bolic fluxes by using the Harten-Lax-van Leer (HLL) pre- 
scription, and damping of the electric field by Ohmic dissi- 
pation that is very stiff by using Strang-splitting techniques; 
Palenzuela et al. (2009, hereafter P09) proposed a numerical 
method that solves hyperbolic fluxes by Local Lax-Friedrichs 
approximate Riemann solver, and the stiff part by using 
implicit-explicit (IMEX) Runge Kutta methods. However, 
these methods use light velocity as the characteristic velocity, 
and their numerical solutions are diffusive when one consid- 
ers problems whose characteristic velocity is much lower than 
light velocity. This indicates that their numerical solutions are 
diffusive in many important high plasma /3 dynamics, and also 
their solutions become highly diffusive when the characteris- 
tic velocity of phenomena is much lower than light velocity. 
In particular, when one solves the dynamics of the accretion 
disk around a black hole with a relativistic jet, one has to 
use relativistic resistive MHD code that can solve both highly 
relativistic and non-relativistic dynamics with resistivity for 
the following three reasons: (1) the saturation of the magne- 
torotational instability (MRI) depends on the resistivity; (2) 
the dynamics of an accretion disk are not ordinarily relativis- 
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tic, especially, the dynamics of the MRJ is sub-Alfvenic; (3) 
the dynamics of the jet are highly relativistic. For these rea- 
sons, previous schemes are diffusive in such phenomena, and 
we need more accurate numerical schemes. We are devel- 
oping a new numerical scheme capable of accurately solv- 
ing problems whose characteristic velocity is quite different 
from light velocity. In this scheme, we obtain numerical flux 
of fluid by using sound velocity as the characteristic veloc- 
ity, and numerical flux of electromagnetic field by using ap- 
propriate characteristic velocities of RMHD. This enables us 
to obtain accurate numerical results when we consider prob- 
lems whose characteristic velocity is much lower than light 
velocity. In addition, P09 pointed out that the Strang-splitting 
method used in the Komissarov method is unstable when ap- 
plied to discontinuous flows with large conductivities. How- 
ever, we find that this problem is not related to the Strang- 
splitting method, but the evolution of electric field E during 
the primitive recovery, that is introduced in the method by 
P09. By considering this procedure, we can apply the Strang- 
splitting method to discontinuous flows with large conductiv- 
ities. 

This paper is organized as follows. In Section |2] the basic 
equations of resistive RMHD are presented. In Section [3] we 
present the numerical method. Results of numerical test prob- 
lems previously presented are shown in Section 2] In Section 
|5] we present results of numerical test problems that cannot 
be solved accurately by previous codes. 



2. BASIC EQUATIONS 
Throughout this paper, we use the units 

c=l, 



(1) 



In Cartesian coordinates, the Minkowski metric tensor ?7^,y is 
given by 

ry,,, = diag(-l, 1,1,1). (2) 

Variables indicated by Greek letters take values from to 3, 
and those indicated by Roman letters take values from 1 to 3. 

2.1. The Maxwell equations 
The covariant Maxwell equations can be written as 

d^Fi"' = Ii', (3) 
d;F'"' = Q, (4) 

where F**" is the Maxwell tensor, *F'^'' the Faraday tensor, 
and the four- vector of electric current. 

If we consider highly ionized plasma, the electric and mag- 
netic susceptibilities can be neglected. Then, one has 



where 



J — 2 ^ pCTi 



Hupa _ r— 



(5) 



(6) 



is the Levi-Civita alternating tensor of space-time, and e^vpa 
is the four-dimensional Levi-Civita symbol. 

We introduce a future-directed unit timelike vector nor- 
mal to a spacelike hypersurface S. Using we can decom- 
pose the Maxwell tensor into following forms: 



(7) 



Similarly, the current four-vector can be decomposed into: 

I>' = qn>' + J>', (8) 



where q is the charge density observed in the rest frame of 
n^, and 7^ the conduction current satisfying J^n^ = 0. In the 
following, we consider only Minkowski spacetime, so = 
(1,0,0,0). 

By using the decomposition of the Maxwell tensor Eq. (jTjl 
and the current four-vector (|8}, the Maxwell equations can be 
spUt into the familiar set 



VB = 0, 
5,E-V xB = -J, 
5,B + V xE = 0. 



(9) 
(10) 

(11) 
(12) 



From Maxwell equations, we can derive the electric charge 
conservation law 



a,^+v-j=o. 



(13) 



2.2. The hydrodynamic equations 

The relativistic hydrodynamic equations can be obtained 
from the conservation of mass, momentum, and energy 

dpN^ = 0, (14) 
d^T^"' = 0, (15) 

where A^^ is the mass density current and T^'' the energy- 
momentum tensor defined respectively as 



where 



N>' = pu'', 

^ "-'fluid"*"-' EM' 



T.Z^phu^^u'^+prj^'', 

TZ = F^^PF;:-h,FP-Fp,)r^>'\ 



(16) 
(17) 



(18) 
(19) 



Here h= 1 + e + p/pis the specific enthalpy, p is the proper 
rest mass density, p is the thermodynamic pressure, and e is 
the specific internal energy. 
The evolution equation of a relativistic resistive MHD is 



' F' ' 



d,{m- j+dj\^F;j I =0, 



(20) 



where D, m', e is the density, momentum density, total energy 
density. In the laboratory frame, D, m, e are given by 



m=p/i7^v + E X B, 
e = ph^^-p+^(E^ + B^), 



(21) 
(22) 

(23) 



where v is the fluid three-velocity, 7 = (1-v^) '/^ is the 
Lorentz factor, and numerical fluxes are 



(24) 



F^j = m'vj + pri'J-E'Ej-B'Bj+-(E^ + B^)r]'-', (25) 



Fl = m'. 



(26) 



This is the most common form of perfect fluid equations for 
the numerical hydrodynamics. 
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2.3. Ohm's law 

The system of Eqs. (|9]l - (fT2] l. ( |20] | is closed by means of 
Ohm's law. Although there are various forms of Ohm's law, 
we consider only the simplest kind of relativistic Ohm's law 
that accounts only for the plasma resistivity, and that assumes 
that it is isotropic similar to previous studies K07 and P09. In 
the covariant form, it is given by 

If^ = aF''''u^ + qou'', (27) 

where cr = 1 /?/ is the conductivity, ij is the resistivity, and qo = 
-I^u'^ is the electric charge density as measured in the fluid 
frame. 

As the Maxwell equations and fluid equations, we can de- 
compose Eq. ( |27] | into 3+1 form, and then the space compo- 
nent of Eq. ( |27] | is given by 

J = cr7[E + vxB-(E-v)v]+^v, (28) 

In the fluid rest frame, Eq. ( |28t becomes 

J = crE. (29) 

The ideal MHD limit of Ohm's law can be obtained in the 
limit of infinite conductivity (cr oo). In this limit, Eq. ( l28T l 
reduces to 

E + vxB-(E-v)v=0. (30) 

Splitting this equation into the components that are normal 
and parallel to the velocity vector, it becomes 

Ei + vxB = 0, (31) 

E||-(E.v)v = 0, (32) 

From these equations, we can obtain the usual result 

E = -vxB. (33) 

3. NUMERICAL METHOD 

In this section, we present our new numerical scheme for 
the resistive RMHD. Since the pioneering studies of resis- 
tive RMHD K07 and P09 use light velocity as the charac- 
teristic velocity, their solution becomes highly diffusive when 
characteristic velocity is much lower than light velocity. In 
our new scheme, we obtain numerical flux of fluid by using 
sound velocity as the characteristic velocity, and numerical 
flux of the electromagnetic field by using Alfven velocity as 
the characteristic velocity. This enables us to obtain accurate 
numerical results even when characteristic velocity is much 
lower than light velocity. In the following sections, we con- 
sider the one-dimensional case. The extension to the multi- 
dimensional scheme using the constrained transpo rt method 
dEvans&HawlevI [19881 Tstone & Norman] [T992I) . wiU be 
shown in our next paper 

3.1. Strang Splitting method 

The relativistic resistive MHD is hyperbolic-relaxation 
equations. In previous work K07 and P09, they assume 
that characteristic velocity is the speed of light. Thus, their 
schemes are highly diffusive when the characteristic veloc- 
ity is lower than light vel ocity. For thi s reason, we apply 
the Strang splitting method ^Strang "1968') and solve the basic 
equations by using each appropriate characteristic velocity. 

First, we split fluid equations Eq. (l20l i as follows: 



where 



d, 




(34) 



K,fiuid = m'v'' + pi\ 



e.fluid 



' m.EM " 



-E'E^-B'B'' + 



^V« = (ExBr. 



(35) 
(36) 
(37) 

(38) 
(39) 



The flux of the fluid component /y/,„y can be calculated by us- 
ing the Riemann solver; the flux of the electromagnetic com- 
ponent F^i^ can be calculated by the method of characteristics. 

Next, we consider the Maxwell equations Eqs. (|9]l - (fT2] i. 
Eqs. ^ and ( fTOl i are not evolution equations but constraint 
equations, and we treat them separately from evolution equa- 
tions. The evolution equations of E and B are Eqs. (fTTT i and 
(fT2l l. By using Ohm's law Eq. (|28] |. Eq. (fTTT i reduces to 



9,E - V X B = -aj[E + V X B- (E • v)v] - q\. (40) 

The source term of this equation includes evolving variables 
E, so this equation is a hyperbolic equation with stiff relax- 
ation terms that requires special care to capture the dynamics 
in a stable and accurate manner Thus, we split the charge 
current J into two parts similar to K07 



Jc = cr7[E + vxB-(E •¥)¥]. 

Then, we split Eq. ( |40] | into two parts 

a,E-V xB = -^v, 
d,E = -J,. 



(41) 
(42) 

(43) 
(44) 



Eq. (03]) is non-stiff equations, and Eq. ( |44T ) is stiff equations. 

As a result, the evolution part of the Maxwell equations can 
be rewritten as 

a,B + VxE = 0, (45) 
5,E-VxB = -^v, (46) 
d,E = -h. (47) 

In component form, Eqs. (|45] | and (|46] | reduce to 

d,B' = 0, (48) 

d,B>'-d,E^ = 0, (49) 

d,B' + d,E>' = 0, (50) 

d,E-' = -qv\ (51) 

d,Ey + d^B^ = -qvy, (52) 

d,E^-d,B^' = -qv\ (53) 

We solve Eqs. ( |49l ), dSOb . ( l52b . and ( |53] | using method of 
characteristics (MOC), which will be shown in Sec. 13.21 Eq. 
( BTl i is solved using the Runge-Kutta method. The numerical 
scheme for the stiff equation Eq. (|47] | will be shown in Sec. 
[331 



3.2. Method of characteristics 

The method of characteristics can be used to solve the initial 
value problems of advective and hyperbolic equations. As is 
well known, the Maxwell equations are hyperbolic, so we can 
solve the Maxwell equations accurately by using this method. 

The Maxwell equations for the transverse fields are Eqs. 
( |49] l, dSOi ). ( |52] |. and ( l53b . By adding and subtracting these 
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equations, for B^, and J-', we obtain 



1 



(54) 
(55) 



where Cd, is the characteristic velocity, and this is equal to the 
speed of light in ordinal Maxwell equations. 
The transverse fields are recovered from by 



B' = ^F--F, 



(56) 
(57) 



The left-hand side of Eq. (l54l i is the total derivative dF/dt for 
an observer moving at velocity ±Cc/,. 

Let us consider conservative discretizations of Eqs. 
and(l52]): 



B"7=B':: + 



At" 

At" r 
Axi 



(py N«+i/2_ 



(fiii/2r'/'-(fiti/2r'/' 



«+l/2 



(58) 



(59) 



where superscript n means the time-step, and subscript ; 
means the coordinate of cell center Using Eqs. (|56] | and (l57t . 
we can obtain the numerical flux of Eqs. ( ISSl l and ( |59] l. (See 
Fig- [I])- The same procedure can be done for time advance of 
E~, B>\ 

The characteristic velocity of the Maxwell equations in vac- 
uum is light velocity. However, since we consider the electro- 
magnetic hydrodynamics equations, appropriate characteris- 
tic velocity has to be used for them. Also, because we con- 
sider resistive systems, the characteristic velocity varies with 
the conductivity a and the scale of wave modes. For example, 
as shown in Appendix, the transverse electromagnetic hydro- 
dynamic waves propagate with the light velocity when k/cr is 
large, where k is the wave number, and they propagate with 
the Alfven velocity when k/a is smaller than a critical value 
depending on p, h, and |Z?| . Because of the finite resistivity, the 
frequency of the transverse waves has an imaginary part cj/ 
(damping rate), which is a increasing function of k/a. In this 
scheme, we use the light velocity as the characteristic velocity 
when cr is smaller than the critical value; when a is larger than 
the critical value, we use appropriate magnetohydrodynamic 
characteristic velocities. The critical value is determined so 
that the transverse waves whose phase velocities are light ve- 
locity are dissipated within the numerical integration timestep 
Af . A detailed procedure to judge whether we use the light 
velocity or magnetohydrodynamic characteristic velocities is 
given in Appendix. 

In addition to the numerical flux of the Maxwell equations, 
the characteristic velocity is also required to construct the 
Maxwell stress terms and the Poynting flux term. When the 
characteristic velocity obtained from the analysis of the trans- 
verse waves is the light velocity, we use the light velocity as 
the characteristic velocity for them; when the transverse wave 
characteristic velocity is the Alfven velocity, we use the char- 
acteristic velocities for them as follows. Note that if the fol- 
lowing characteristic velocities are not used, numerical inte- 
gration becomes unstable, and unstable numerical oscillation 
occurs. 

Then, the necessary procedures are as follows: 




1/2 




n+1/2 n + 1/2 



E V + B7 Ev " B; 



z i + 1/2 



Fig. 1. — A schematic drawing of Eulerian-like chai'acteristics when one 
uses piecewise hnear interpolation, c^h is the characteristic velocity. On the 
left is the subsonic case, and on the light is the supersonic case. These fig- 
ures show that half time-step transverse electromagnetic field E^' and B~ are 
determined by the fields at the base of two characteristics. 



1. For the numerical flux of Maxwell equation Eqs. 
and ( |59] l. we use the Alfven wave velocity in labora- 
tory frame because the information of transverse elec- 
tromagnetic fields are transmitted by the Alfven wave ' 
. In relativistic MHD, the Alfven velo city in labora tory 
frame val can be obtained by solving d Anile II 1 9901) 



-B^ = 0, 



where H = ph + b^, a = ^(val - v''), B = lf- 
is the covariant magnetic field defined as 



b^' = 



B 

7V-B, — + 7(vB)v 

7 



(60) 
VAiy^, and 

(61) 



2. For the Maxwell tension terms -EE-BB in Eq. 
we use the Alfven wave velocity in fluid comoving 
frame vac because magnetic tension force is originated 

' In the case of small limit, the Alfven velocity in laboratory frame 
Val becomes v.,. We find that when this I'v is also small, numerical oscil- 
lations occur and numerical integration becomes occasionally unstable. We 
can prevent this purious oscillations, if we use the characteristic velocity as 
Vf^i when v^/, < 0.1 |-B|/ \/ ph-V |B|-. In this case, our scheme for solving 
the Maxwell equations becomes equivalent to the HLLE scheme. Note that 
introduction of this modification does not change any results presented in this 
paper 
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by the Alfven wave. In relativistic MHD, the Alfven 
velocity in fluid comoving frame is given by 

Vac = —^- (62) 



H 

3. For the Poynting flux E x B of energy equation Eq. ( [39] l 
and the Maxwell pressure terms /2 + /2 in Eq. 
(l38l l. we use the fast magnetosonic wave velocity in lab- 
oratory frame because the magnetic pressure originates 
in the magnetosonic wave. In relativistic MHD, fast 
magnetosonic wave velocity in the laboratory frame Vf,„ 
can be obtained by solving 

ph(l-c^Ja^ = (l-v%,)m' + phc'>^-cyi (63) 
Eq. ( l63T l is a quartic equation, and in an ordinary one 
has to use the Newton-Raphson method or the quar- 
tic formula for obtaining solutions. However, since 
our scheme splits the fluid part and the electromagnetic 
part, the sound velocity c, can be set equal to zero. 
Then, the characteristic equation Eq. (l63T l reduces to 

ph^Hv'-Vf,„f = (l-v}j\b\\ (64) 
By using the quadratic formula, one can obtain solu- 
tions of above equation: 

p/i7V ± \b\ ^y\b\^ + (l- 



ph"f^+ \b\'- 



(65) 



To sum up, we only have to substitute the appropriate char- 
acteristic velocities val,vac, and vy,„ into q/, in Eq. ( |54] l, 
and calculate the electromagnetic field E,B at half time step. 
Then, the numerical fluxes of electromagnetic hydrodynamics 
equations are given by 

'1 



^m.EM - '^Ac'^Ac 



-B'acB\c- 



(66) 

/, 

X Bf„,y\ (67) 
where Eac,Bac means that they are calculated by using the 
Alfven velocity in comoving frame, and Ef,„,Bf,„ by using 
the fast magnetosonic wave velocity in laboratory frame. For 
the numerical flux of the Maxwell equation, one has to use the 
Alfven velocity in laboratory frame val for the calculation. 

3.3. Stijfpart 

As explained Sec. 13.11 Eq. (l44l i contains stiff terms. Fol- 
lowing the previous work K07, we split the equation into com- 
ponents normal and parallel to the velocity vector. 

9,E||+cr7[E||-(E-v)v] =0, (68) 
9,E_L + cr7[E_L + vxB)]=0, (69) 
Since we use the Strang splitting method, the right-hand side 
of the above equations can be considered constant other than 
the electric field E. As a result, these equations can be solved 
analytically 



E 



= E|j exp 



a 

— t 

1 . 



(70) 



Ei = El + (E'i-El)exp[-a7f], (71) 
where E*^ = -v x B and suffix indicates the initial compo- 
nent. If we use the explicit integrator, the stiff equation has to 
be solved in very small time steps Af. However, since Eqs. 
( ITOI i and dTTI ) are formal solutions, we can avoid the stability 
constraints of the time step. In the context of ambipolar diffu- 
sion in partially ionized plasma, a similar numerical technique 
using the piecew ise formal solution of stiff part is known to be 
usefu l scheme (llnoue et al. 1 120071; llnoue & Inutsukall2008l 
l2009l) . 



3.4. Constraint Equations 

It is well known that Eqs. Q and ( fTOl l are constraints on the 
Cauchy surface. Though Maxwell equations ensure that these 
constraints are preserved at all times, straightforward numer- 
ical integration of Maxwell equations does not preserve these 
properties because of the accumulated numerical error This 
causes corruption of numerical results, and results in a crash 
in the end. For this reason, there are a number of numerical 
techniques for avoiding this problem. We have implemented 
hyperbolic divergence cleaning for the electric field. The main 
idea of the hyperbolic divergence cleaning is that one defines 
new variable ^I^ as the deviation from constraint equations, 
and arranges a system of equations to decay or carry the de- 
viation out of the computational domain by high speed 
waves. For the magnetic field, if one sets constant, the 
constraint equation can be satisfied in one-dimensional case. 
In the multi-dime nsional case, we can imp l ement constrained 
transp ort method (Evans & Hawley I I1988I: IStone & Norman I 
|1992b . The detailed implementation will be presented in our 
next paper 

For hyperbolic divergence cleaning, we modify Eqs. (|9} 
and dsTl ) 



d,E'' + d,^ = -qv\ 



(72) 
(73) 



where is a new dynamic variable and k a positive constant. 
Clearly, when we set = 0, we can recover standard Maxwell 
equation Eq. (|9|l. From these equations, we can obtain the 
telegrapher equation for 



a?*+Kai'-v^* = o. 



(74) 



Thus, propagates at the speed of light and decays exponen- 
tially over a timescale 1 /k. 

Similar to Eq. ( |44] |. Eq. (l72T i contains stiff source terms. 
Thus, we split the equation into a stiff part and non-stiff part 



a,^' + V-E = ^, 

d,'^=-K^. 

The analytical solution of Eq. ( 176] ) is 

* = ^'oexp[-Atf] 

where is the initial value of 

3.5. Primitive recovery 



(75) 
(76) 

(77) 



In order to compute numerical flux 
and ( [39] l, the primitive variables {p,v,/?,B,E} have to be re- 
covered from the conserved variables {D, m, e, B, E}. In con- 
served variables, E and B can be obtained by evolving the 
Maxwell equations. However, as pointed out by P09, it is 
more stable to perform evolution of stiff part Eqs. (iTOt and 
( fTll i during this primitive recovery process when a is large, 
i.e., ideal MHD approximation is valid. This is because when 
we consider MHD approximation, the electric field E is equal 
to -V X B; however, in general, primitive recovered E does not 
satisfy this relation. In what follows we explain the primitive 
recovery procedure following P09. 

1 . Set an initial guess for the velocity by using previous 
time step value Then, evolve electric field E using Eqs. 
dTOll and dzB- 

2. Subtract Poynting flux and electromagnetic energy den- 
sity from conserved variables, and new variables can be 
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defined as follows: 

m' = ph-fh, (78) 

e' = ph-f^-p. (79) 

Then, variables {D,m',e'} are the ideal relativistic fluid 
conserved variables, and can be recovered by using the 
ordinary procedures. 

3. Replace the initial guess for the velocity with the ob- 
tained velocity v, and repeat the steps 1-3 until the 
primitive variables converge. 

3.6. Algorithm 

In this section, we provide the detailed numerical algo- 
rithm. 

In Cartesian coordinates, the relativistic resistive MHD 
equations written in conservative fashion are 




where 



m = phj^\ + E X B, 



- e, fluid ■ 



, = -E'E'-B'B' + 



(80) 

(81) 
(82) 

(83) 

(84) 
(85) 
(86) 

(87) 
(88) 



^;% = (ExB)\ 

The electric field E and magnetic field B are evolved by the 
Maxwell equations. If the Ohmic dissipation is considered, 
the Maxwell equations have stiff and non-stiff part. The non- 
stiff part is 

5rt/Maxwell + ^.c^Maxwell = '^non-stiff , (89) 



Maxwell ~ 



^non— stiff ~ 



B> 




( ^ \ 


B" 




Ey 


E' 
E>- 


J ^Maxwell — 


B~ 


E~ 




-By 






E' 






\r ) 


( ° 




\ 















-qvy 








q 

\ 


) 





(90) 



(91) 



where F = alE" + (v x B)-" - (E • y)v'] + qv\ The Maxwell 
equations are consistent with the equation of charge conserva- 
tion. However, numerical errors in general destroy the conser- 
vation law in a way similar to the constraint equations. Thus, 



the above equation contains the equation of charge conserva- 
tion. 

As explained in Sec. I3.3l and [34l the stiff -part is evolved by 
using the formal solution 



E||=E|I exp 



a 

— t 

7 . 



(92) 



Ei=El+(E'i-El)exp[-a7f], (93) 

* = *oexp[-Kf]. (94) 

Using the above system equations, the second-order numer- 
ical algorithm is given as follows. 

1. Advance the Stiff -part equations over Af/4 by using 
the formal solutions Eqs. (|92] | - (|94] |. 

2. Advance the non-stiff part of Maxwell equations Eqs. 
(|89] | and W\\ over Af/2 by using method of char- 
acteristics as explained in Sec. 13.21 and calcu- 
late numerical flux Fem ( 1381 ) and ( [39] l. On the 
other hand, numerical flux Ffiuid ( |35] | - dJTl l can 



jMarti&MUner 1994; Marti 1996; Banvuls et al. 


1997 


; Alov et al. 1999; Pons et al. 2000; Font et al. 


200C 


; |Del Zanna & Bucciantini 2002|: Marti 


MiiUer 


2003 


; iMianone & Bodo n2005i; iMianone et al. 


12005b. 



In this paper, we use the HLLC solver. 

3. Advance conserved variables D, m, e over the half time- 
step Af/2 by using Eqs. ( |34] | - ( [39] l. Then, calculate 
primitive variables of half time step {/""""'/^ by primitive 
recovery explained in Sec. 13.51 In our scheme, electric 
field E has to be evolved Af /4 by using formal solution 
(|92| i and (|93] l during primitive recovery. Primitive vari- 
ables obtained through this procedure are used for the 
calculation of the numerical flux at f = f + Ar/2. 

4. Again, advance initial stiff variables over Af/2 by us- 
ing the formal solutions Eqs. ( |92] | - (|94] |. 

5. Calculate temporal second-order numerical flux (l35i - 
(|39] | by using primitive variables obtained through the 
procedure 3. Then, advance conserved variables D, m, e 
over Af by Eq. ( l34l l. and electric field E and magnetic 
field B by the Maxwell equations of stiff-part Eq. dgTT l. 

6. Calculate primitive variables by a primitive recovery 
process. During this process, the electric field E is ad- 
vanced over Af/2 by using formal solutions (l92l l and 
(193). 

For the spatial second-order, we use the MUSCL scheme 
by Van Leer explained in Appendix. 

Note that if we evolve electric field E in integration of stiff 
equations or primitive recovery procedure, we have to evolve 
other primitive variables. This is because conserved variables 
are not changed during those procedures, and this means that 
the change of electric field E affects all other primitive vari- 
ables. 

4. TEST SIMULATIONS 

In this section, several one-dimensional test simula- 
tions given in previous studies K07 and P09 are pre- 
sented. For the numerical flux of fluid, we use the HLLC 
solver (iMignone & Bodo1l2005h . We use an ideal equation 
of state pe = p/iT- 1) with F = 2, and Courant number, 
CFL = 0.25. 
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4. 1 . Large amplitude CP Alfven waves 

This test consists of the propagation of a large ampli- 
tude circularly -polarized Alfven waves along a uniform back- 
ground field Bo- The analytical exact solution of this problem 
is given by Del Zanna et al. (2007) jDel Zanna et al.ll2007h : 
and this problem is used as the ideal-MHD limit test problem 
by P09. We use the same condition as P09. 

{By,B') = luBo (cos[fe(x- VAf)], sin[^(x- VAf)]) , (95) 
(v^',/) = -^(B>,B^), (96) 

where = Bq, v' = 0, A: is the wave number, and 77^ is the 
amplitude of the wave. The special relativistic Alfven speed 
va is given by 

For the initial data parameters, we have used p = p = r]A = I, 
and B() = 1.1547. Using these parameters, the Alfven velocity 
is V/i = 1/2. For the boundary condition, the periodic one is 
used. In addition, we use a high uniform conductivity a = 
10* following P09, since this is the exact solution of ideal 
relativistic MHD. 

Fig. 12] is results of our new code at f = 2.0 (one Alfven 
crossing time) for three different resolution cases with = 
{50, 100,200}. The computational domain is x e [-0.5,0.5]. 
This result indicates that our new code reproduces ideal rela- 
tivistic MHD solutions when the conductivity a is high. 
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Fig. 3 . — Li norm eiTors of the tangential magnetic field By under different 
grid resolution for the second-order schemes using the new scheme. The 
left-hand side is the result of Large amplitude CP Alfven waves, and the 
light-hand side is the result of the self-similar current sheet. 
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Fig. 2. — The results of large amplitude circularly-polaiized Alfven wave 
test with large conductivity a = 10*. This test is carried out for three different 
grid points A' = 50, 100.200. 

In these test problems, we cannot achieve full second-order 
accuracy. The left-hand side of Fig. |3]is the Li norm errors 
of the tangential magnetic field By of this test problem. This 
figure shows that our numerical result is nearly 1 .5-order con- 
vergence. We estimate this is because our scheme uses many 
operator splittings, and the time accuracy of our scheme wors- 
ens. Note that this problem is one of the most difficult to solve 
in relativistic resistive MHD, since this is the limit of large 
conductivity cr. The right-hand side of Fig. [3]is the Li norm 
errors of the B, of the next test problem. Since the conductiv- 
ity (T is moderate value in that test problem, our new scheme 
achieves second-order convergence. 

4.2. Self-similar current sheet 



This problem is used as the test problem of highly resis- 
tive cases in K07 and P09. In this test, it is assumed that the 
magnetic pressure is much smaller than gas pressure, so that 
the background fluid is not influenced by the evolution of the 
magnetic field. We assume the magnetic field has only tan- 
gential component B = (0,B(x,f),0), and B(x,t) changes its 
sign within this current sheet. Since we are interested only in 
the evolution of the magnetic field, the background fluid is set 
initially in equilibrium, p = const. In addition, we assume that 
the conductivity a is high, and the diffusion timescale is much 
longer than the light propagating timescale. Although the re- 
sistive relativistic MHD equation is hyperbolic, this assump- 
tion allows us to neglect the displacement currents at least in 
the rest frame. As the result, the evolution equation is reduced 
to 

d,B--d^B = 0. (98) 

(7 

This equation has exact solution 

B(x,f) = BoerfQy|^, (99) 

where erf is the error function. Following K07 and P09, we 
set the initial condition at f = 1 with p = 50, p= 1, E = v = 0, 
and cr = 100. The computational domain is [-1.5, 1.5], and 
the number of grid points is N = 200. Fig. |4]is the numerical 
result at t = 9. This figure shows that our scheme can solve 
a highly resistive problem accurately. The convergence rate 
is consistent with the second-order spatial and temporal dis- 
cretization. 
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-1.5 -1 -0.5 0.5 1 1.5 



X 

Fig. 4. — The result of self-similar current sheet test comparing the ex- 
act solution. The solid line is the exact solution, and the dotted line is the 
numerical result with conductivity cr = 10". 

4.3. The propagation of Aljvenic transverse waves with 
Ohmic dissipation 

In order to confirm the capability of our method for the rel- 
ativistic resistive MHD, we perform the test calculation of the 
propagation of Alfvenic transverse waves with Ohmic dissi- 
pation, and compare the results with the exact dispersion re- 
lation Eq. ( lAlOb obtained in Appendix. 

As explained in Appendix, the resistive relativistic mag- 
netohydrodynamic equation contains transverse wave modes 
that become the light wave in large k/a region, and become 
the Alfven wave in small k/a region. To demonstrate the 
propagation of transverse waves, we set the initial condition 
by e igenfunctions of the mode obtained from Eqs. ( IA6b - 

B' = 0.05cos(kx), (101) 




£^ = ^fi^ (104) 

where w' = w/cr and k' = k/ a = 2tt / a, and uj is the solution of 
the dispersion relation Eq. (lAlOl l. We set the same parameters 
in Appendix. 

(p/i,fi^) = (1.5,0.55). (105) 

Since the enthalpy includes the information of the equation of 
state, one can take any value of F. In this calculation, we set 
r = 2 and p=\. The computational domain covers the region 
[-0.5,0.5] where the periodic boundary condition is imposed, 
and the number of grid points is N = 200. 

The propagation speed of the numerical waves can be de- 
termined by tracing the position where B, is maximum. We 
measure the propagation speed and evaluate Re[w] based on 
the time when the maximum of reaches x = again, i.e. 
one-wave crossing period. The damping rate Im[cj] is mea- 
sured by using Z?]^ = Z?j^)exp[Im[w]f] where is the maxi- 
mum of B" after the one-wave crossing time. 

In Figs. |5] we plot the real and imaginary part of u/a 
against k/a. The solid line is the exact dispersion relation 
obtained in Appendix. We have performed the calculation in 
the cases of /t/a = 0.01,0.1,0.5, 1,4, 10, 100. These figures 
show that our new numerical code can reproduce the propa- 
gation of Alfvenic transverse waves accurately for any value 



of the conductivity a. 
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Fig. 5. — The result of the propagation of Alfvenic transverse waves with 
Ohmic dissipation test problem. The solid line is the exact dispersion relation, 
and the dots ai'e the numerical solutions. The test calculations are performed 
in the cases of k/a = 0.01,0.1,0.5, 1,4, 10, 100. 



4.4. Shock-tube problem 

For the shock tube test problem, we consider the simple 
MHD version of the Brio and Wu test as P09. The initial left 
and right states are given by 

(P^,/,(bY)= (1.0,1.0,0.5) for jc<0.5 (106) 
(/,p*,(B-'/) = (0.125, 0.1, -0.5)for x > 0.5 (107) 

All the other fields are set to 0. 

Fig. |6]is the numerical results at f = 0.4 that change grid 
points = 100,200,400. The computational domain covers 
the region [0, 1]. We also plot an ideal RMHD solution by 
the solid line computed by a publi cly available code devel- 
oped by Giacomazzo and Rezzolla dGiacomazzo & RezzoUa I 
|2006|) . The conductivity is uniform with a = 10^. The solu- 
tion of this Riemann problem contains a rarefaction moving 
to the left, a shock moving to the right, and a tangential dis- 
continuity between them. Fig. |6] shows that our numerical 
solution of the resistive MHD can reproduce the profile of an 
ideal MHD shock tube problem using high conductivity cr. In 
addition, our numerical solution captures contact discontinu- 
ity as sharp as P09. 

Fig. |7] is the numerical results of the same problem that 
changes the conductivity cr = 0, 10, 10^, 10^, 10^. We also plot 
the ideal RMHD solution by the solid line. The number of 
grid points is N = 400. This result shows that our numerical 
solution reproduces nearly the same results as P09. 

P09 reports that Strang's splitting method becomes unstable 
for moderately high values of the conductivity for this shock 
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tube problem, and one has to use the implicit method. How- 
ever, this is not related to whether one uses Strang's splitting 
or implicit method, but to the revision of the electric field dur- 
ing the iteration of the primitive recovery (H. R. Takahashi 
2010, private communication). Our scheme uses Strang's 
splitting, but can solve this shock tube problem stably even 
when a > 10^ if we revise the electric field during the primi- 
tive recovery as explained in Sec. 
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Fig. 6. — The numerical results of the Riemann shock tube test problem 
for three different grid points N = 100,200,400. We use the conductivity 
a = 10*. The sohd Une is the ideal solution. 



0.6 



0.4 

0.2 


-0.2 
-0.4 
-0.6 
-0.8 



ideal 
c = 
a= 10 
a = 10, 

o = iaJ 

= 10^ 



0.2 



0.4 



0.6 



0.8 



Fig. 7. — The numerical results of the Riemann shock tube test problem 
for different conductivity cases: o- = 0, 10, 10-, 10^, 10'^. The number of grid 
points is N = 400. The solid line is the ideal solution. 



5. TEST SIMULATIONS FOR FLUID DOMINATED 
CASE 

The previous studies K07 and P09 use light velocity for the 
characteristic velocity. Thus, their numerical solutions be- 
come highly diffusive when one considers problems whose 
sound velocity or Alfven velocity is much lower than light 
velocity. In this section, we perform test problems in such 
cases, and compare the results of the HLL code with that of 
our code. 

5.1. Shock tube test problem 

In this section, we compute a high plasma /3 shock tube 
problem, and compare results of our code with those of the 
HLL code. The initial left and right states are given by 

(/,//,(B0^)= (104,1.0,0.05) for jc<0.5, (108) 
(/,/,(£-'/) =(10^0.1,-0.05)for x>0.5. (109) 



All the other fields are set to 0. 

Figs. [8] are the numerical results of our code and the HLL 
one, being compared with ideal solutions at f = 30.0. The 
number of grid points is N = 400. These figures show that 
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Fig. 8. — The numerical result of shock tube problem for fluid energy 
dominated case comparing with that of HLL and ideal solution. On the left 
is the density profile, and on the right is the profile of the tangential magnetic 
field . The number of grid points is N = 400. 

HLL solver becomes more diffusive than our code. In addi- 
tion. Figs. [8]show that the density profile of the shock heated 
region somewhat overshoots that of the ideal solution, and 
tangential magnetic field slightly undershoots that of the 
ideal solution. These results show that when the plasma f3 is 
high, the HLL solver becomes highly diffusive and does not 
reproduce the correct value of the shock heated region. In 
contrast, our numerical results reproduce ideal solutions very 
well even for high f3 problems. 

5.2. The propagation of contact discontinuity 

In this section, we calculate the propagation of a contact 
discontinuity, and study the accuracy of capturing the con- 
tact discontinuity for various advection velocities. When one 
uses the HLL code by Komissarov, the numerical results can 
be expected to be diffusive for the case of very slow advec- 
tion velocity, since the HLL code uses light velocity for the 
characteristic velocity. In contrast, our new code uses sound 
velocity for the fluid characteristic velocity, and the numerical 
results will be more accurate for any advection velocity. 

We consider the propagation of contact discontinuity of 
magnetohydrodynamics. The initial condition is 

(p^,/,(BO^)= (1.0,1.0,0.1) for x<0, (110) 
(p*,/?*,(B-'/) = (1.5,1.0,0.05)for x > 0. (Ill) 

All the other fields are set to 0. 
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Since we want to consider the ideal fluid case, we con- 
sider high conductivity cr = 10^. We use an equation of state 
with r = 5/3, and the computational domain covers the region 
[-0.5,0.5] with 100 grid points. The CFL number is 0.25, and 
the integration is carried out until 2 fluid crossing time. For 
the boundary condition, the periodic one is used. For the ad- 
vection velocity, we use the following velocities: 



v'=0.9, 0.5, 0.1, 0.05, 0.01. 



(112) 
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Fig. 9. — The numerical results of the propagation of contact discontinuity 
of RMHD by using our new scheme for different advection velocity: v' = 
0.9,0.5,0.1, 0.05, 0.01 . The number of grid points is N = 100. 
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Fig. 10. — The numerical results of the propagation of contact discon- 
tinuity of RMHD by using the HLL code for different advection velocity: 
= 0.9, 0.5, 0. 1 , 0.05, 0.01 . The number of grid points is N = 100. 



characteristic velocity. The diffusive result of HLL code is al- 
ways problematic for any discontinuity when the propagation 
velocity is much smaller than light velocity. In contrast, since 
our code uses appropriate characteristic velocities, the numer- 
ical dissipation does not depend on the characteristic velocity. 
For this reason, our new code can solve any advection veloc- 
ity problems accurately, especially for the problems including 
discontinuities. 



On the left of Figs. |9]are the numerical results of the den- 
sity profile calculated by using our new code, and on the left 
of Figs. [TO]are the numerical results of the density profile cal- 
culated by using the HLL code. The solid lines are the ideal 
solution. This figure shows that the numerical results of den- 
sity profiles by HLL code of = 0.9 is nearly equal to that of 
our new code. However, the numerical results by HLL code 
become more diffusive than by our code as the advection ve- 
locity becomes small; in contrast, the accuracy of the numeri- 
cal results by our code is nearly independent of the advection 
velocity. The right hand side of Figs. |9]are the numerical re- 
sults of tangential magnetic field B^' calculated by using our 
new code, and the right hand side of Figs. [TO]are the numeri- 
cal results of tangential magnetic field B'' calculated by using 
the HLL code. Similar to the density profile, the numerical 
results by using the HLL code become more diffusive than by 
using our code as the advection velocity becomes small; in 
contrast, the accuracy of the numerical results by our code is 
nearly independent of the advection velocity. 

In conclusion, the HLL code is not capable of accurately 
solving problems whose advection velocity is smaller than 
Ught velocity, since the HLL code uses light velocity for the 



5.3. The propagation of small amplitude Alfven wave 

In this section, we consider the propagation of small am- 
plitude Alfven waves in high /3 plasma. The integration is 
performed for different resolutions, and we compare the nu- 
merical results of Komissarov's HLL code and our code. For 
the application to the numerical simulation of MRl, the in- 
tegration is performed for a small number of grid points: 

= 16,32,64 for one wavelength of the Alfven wave; this 
corresponds to the number of grid points for resolving the 
wavelength of maximum growth rate of MRl. 

For the initial condition, we consider 

(p,/?,B-\B-'') = (10,0.05,0.1,0.1), (113) 
B-' = 0.01sin(27rjc/L), (114) 



In this case, the Alfven velocity va and plasma beta /3 are 
given by 

va = 3.14x10"^ /3 = 5.02xlO^ (116) 
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where the Alfven velocity and the plasma beta is defined as 



VA = 



' \B\ 



(117) 
(118) 



Since the initial magnetic field is very weak for most of the 
MRl phenomenon, a weak magnetic field is considered. In 
order to consider the ideal fluid case, we set a high conduc- 
tivity (7 = 10*. We use an equation of state with F = 2. The 
computational domain covers the region [-0.5,0.5]. The CFL 
number is 0. 1, and the integration is carried out until 1 Alfven 
wave crossing time. For the boundary condition, the periodic 
one is used. 
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Fig. 1 1 . — The numerical result of the propagation of a small slow Alfven 
wave. On the left is the result of our new scheme, and on the right is that of 
the HLL code. The number of grid points is A' = 16, 32, 64. 



The numerical results of our code and HLL are presented in 
Figs. [TT| Although the amplitude of both results falls because 
of the numerical diffusion, it can be seen that HLL results are 
more diffusive than our numerical results when the number of 
grid points is N = 16, 32. When the number of grid points is 

= 64, the numerical result of HLL code is a little more ac- 
curate than that of our code. This is because our new scheme 
uses an operator split for the accuracy, and the convergence 
rate is a little less than second order in time. However, from a 
practical point of view, it is impossible to cost 64 grid points 
for the wavelength of maximum growth rate of MRl in many 
cases, and still our new code can integrate the growth of mag- 
netic field by MRl more accurately. 

In conclusion, when one considers the high /3 plasma, our 
code is more accurate than the HLL code because our code 
uses sound velocity and Alfven velocity as the characteristic 
velocity. In particular, the above results show that our new 



method is useful for the application to the phenomena includ- 
ing MRl. This instability occurs in the system whose angular 
momentum changes as r""(0 < n < 2), and the amplitude of 
the perturbative Alfven waves grows exponentially in over the 
duration of nearly one Kepler rotation. Since the above condi- 
tion is satisfied in most of the differential rotating systems in 
gravity, MRl is one of the most important astrophysical phe- 
nomena. In order to reproduce this instability numerically, 
one has to resolve the wavelength of maximum growth rate. 
However, this is difficult for most problems, since this wave- 
length is proportional to the initial weak magnetic field. For 
this reason, in order to reproduce MRl numerically, one has 
to use numerical schemes that can integrate small amplitude 
Alfven waves accurately by smaller number of grid points. 
Then, the results of test problems in this section show that 
our new numerical scheme can deal such problems more ac- 
curately than previous codes. 

In these three test problems, we consider extremely high 
density cases in order to distinguish differences easily. How- 
ever, this can always happen when the magnetic field is weak. 
As a result, if one considers problems including an initially 
weak magnetic field like magnetic rotational instability (MRl) 
in the accretion disk, our code can produce more accurate re- 
sults. 

6. CONCLUSION 

In this paper, we have presented a new numerical scheme 
of resistive RMHD for one-dimensional case which can solve 
matter dominated problems more accurately than the exist- 
ing numerical method. Since this new scheme uses differ- 
ent characteristic velocity for obtaining the numerical flux of 
fluid and electromagnetic field, one can solve accurately and 
stably problems whose characteristic velocity is much lower 
than that of light. 

When one considers relativistic problems, one has to solve 
stiff equations for electric fields. In general, it is difficult to 
deal with stiff equations, and special methods have been pre- 
sented; for example, K07 uses the Strang's splitting method, 
and P09 use the implicit method. P09 report that Strang's 
splitting method is incapable of solving problems that include 
discontinuity, such as shock. However, we find that this is not 
related to the method for the stiff equations, and one can solve 
problems including shock if one evolves the electric field dur- 
ing the primitive recovery; we use Strang's splitting method, 
and the solver is well behaved for shock tube problems. The 
results of other test problems show that our new scheme is ca- 
pable of accurately solving both highly resistive problems and 
nearly ideal MHD ones. In addition, it has been shown that 
our code can solve low characteristic velocity problems more 
accurately than the HLL code. 

The problems of high density and high plasma /3 appear 
when one considers magnetorotational instability (MRl) in 
the accretion disk with a relativistic jet, for example. In this 
case, one has to use relativistic resistive MHD code that can 
solve both highly relativistic and non-relativistic dynamics 
with resistivity for the following three reasons: (1) the satu- 
ration of the magnetorotational instability (MRl) depends on 
the resistivity; (2) the dynamics of an accretion disk are not 
ordinarily relativistic, especially, the dynamics of the MRl is 
sub-Alfvenic; (3) the dynamics of the jet are highly relativis- 
tic. Our new scheme can solve such problems accurately even 
when the initial magnetic field is very weak. 

We present multi-dimensional extension of our scheme in 
our next paper 



12 



TAKAMOTO & INOUE 



We would like to thank Bruno Giacomazzo for provid- 
ing the code computing the exact solution of the Riemann 
problem in ideal MHD. Numerical computations were in part 
carried out on Cray XT4 at Center for Computational As- 



trophysics, CfCA, of National Astronomical Observatory of 
Japan. This work is supported by Grant-in-aids from the Min- 
istry of Education, Culture, Sports, Science, and Technology 
(MEXT) of Japan, No. 22-3369 (T. I.). 



APPENDIX 

THE DISPERSION RELATION OF THE RELATIVISTIC ELECTROMAGNETIC FLUID 

As explained in Sec. 13.21 we solve the evolution of electric and magnetic field by the method of characteristics. For the 
characteristic velocity, we use the appropriate MHD characteristic velocity when a is large, that is, ideal MHD approximation is 
valid. However, when the conductivity a is not so large, we have to replace the characteristic velocity with the speed of light. In 
this section, we discuss when to switch the characteristic velocity from appropriate MHD characteristic velocity to light velocity. 
In the following, we calculate the linear perturbation of the relativistic electromagnetic equation in order to obtain characteristic 
velocity. 

The relativistic electromagnetic fluid equations are given by 

p/iM^a^u'=-V/? + (^E+JxB), (Al) 
a,B = -VxE, (A2) 
5,E = VxB-J, (A3) 
J = cr7[E + vxB-(E-v)v] + 9V, (A4) 
g = V-E, V-B = 0. (A5) 
To obtain the dispersion relation, we start by expanding physical variables around an unperturbed state in the following frame: 

• The fluid is at rest: Vo = 

• The x-coordinate is parallel to the k: k = k^,^ 

• The magnetic field is in the x-direction: Bq = B^e^; 

• charge neutrality: = 0, Eq = 

Since we only want to judge when to switch characteristic velocity, we consider propagation of the transverse waves along the 
magnetic field. When one uses this procedure during the numerical simulation, one only has to calculate B^-E^ of the simulation 
data, and substitute its square root into the above B". This is because B^-E^ is a scalar, and becomes the square of the magnetic 
field in the fluid comoving frame because of the assumption of the charge neutrality. Since the magnetic field appears only in 
the form of B^ in the following procedure, one can neglect the sign of magnetic field. In the following, we consider only the 
characteristic velocity of transverse waves. 

In the above condition, the Alfven mode is included in the z-component of the velocity 5v^ and magnetic field 5B~, and 
decouples from other variables. For this reason, we consider only variables related to 6v~ and 5B^. 

We replace the current vector in Eq. (lAU with Eq. (IA3) . Then, the perturbed equations are 

iujphSv- + ikB'SB' - iujB'SE^' = 0, (A6) 

iujSB'-ikSE^' = 0, (A7) 

aB'Sv' + (iuj - (j)5E' = 0, (A8) 

aB-'Sv' + ikSB' + (a-iuj)SE>' = 0. (A9) 
From these equations, the following dispersion relation is obtained: 

phuj'^ + ia(B^ + 2ph)uj^ - [k^ph + a\B^ + ph)]up- - ia(B^ + ph)k^uj + a^(B')V = 0, (AlO) 

Eq. (lAlOl i is the biquadratic equation with respect to uj, and has the formula of radicals. However, the analytical formula is 
very complex and hard to analyze, and is not suitable for obtaining the characteristic velocity. As shown below, the transverse 
wave becomes an A lfven wave in the long wavelength regime, and the light wave in the short wavelength regime shown in 
Fig. lAll Fig. lAll shows that the damping rate is a monotonically increasing function of k. For this reason, we establish the 
following criterion for the characteristic velocity: when all light modes damp during one time step At" we use appropriate MHD 
characteristic velocity for the method of characteristics; when some light modes do not damp during one time step Af", we use 
light velocity for the method of characteris tics. W e will discuss this method in detail in the following. 

First, we substitute uj = ujR + iuj into Eq. dAlOb . and divide the dispersion relation into a real part and imaginary part. Then, the 
real part is 

- [phk^ + a^(B^ + ph) + 3a(B^ + 2ph)uj,+6phujj]ujl 
+ phujj + a(B^ + 2ph)uj] 

+ [phk^ + a^iB^ + ph)]ujj + aiB^ + ph)k^uJi + iB''fa^k^ = 0, (All) 

and the imaginary part is 

[B^cr + 2p(cr+2tj/)]a;«-[p/z((T+2tj/){/t^ + 2a;/(CT + a;/)} + BV{/t^ + a;/(2(T+3a;/)}]wfi = 0, (A12) 
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Eq. ( IA12b implies that the solution for ljr is and a conjugate complex numbers. The solutions of = are pure decaying 
modes, so the other modes a re th e desired propagating ones that become light velocity in the limit of small a and Alfven velocity 
in the limit of large cr. Figs. lAll are the dispersion relation for the propagation modes of the following parameters: 

(p/i,5-0 = (1.5,0.55) (A13) 

These figures show that this mode becomes light in the limit of small cr and Alfven wave in the limit of large a. Although the 
form of ojn does not become as Fig. lAll for some parameters, this mode always becomes light waves in the limit of small a. 




0.001 



— — light mode 
Alfven mode 



1 

k/a 




0.001 



0.0001 



1e-05 



100 



Fig. A1. — The dispersion relation of tlie propagation mode. The left-hand side is the real part of oj/cr, and the right-hand side is the imaginary part of u/cr. 
In the figure of the real part of w/cr, we also plot reference lines whose phase velocities are Alfven velocity (long-dashed line) and speed of light (short-dashed 
fine). The phase velocities can be obtained from the data al k/a = 0.01 by using the formula Cpi,„„ = ui/k. The parameters are set as (ph.B'') = (1.5,0.55), and 
the Alfven velocity is given by va ~ 0.409. This figure shows that this mode becomes hght waves in the limit of large k/a and Alfven waves in the limit of small 
k/a. In addition, this mode has a maximum decay rate in the limit of large k/a. 



From Eq. ( IA12l i. this desired mode can be obtained as follows: 

2 ph(a + 2Lo,) [k^ + 2w/(cr + uj,)]+ B^a [k^ + w/(2cr + 3w/)] 



B^a + 2ph(a+2uji) 



We substitute this into Eq. ( lAl 11 1, and we obtain 



where 



(A14) 



a4k +a2k +ao = Q, (A15) 



a4 = -B^p^h^aia + lu,) - p\^(a + 2oj,f , (A 1 6) 

a2 = B'^iB'fa* -B^a\a + 2uj,)- 2B^pha^(2c7^ + 1(juj, + 6ujj) 

- p^h\a + 2uj,f{-4-(B'fa^ + 2ph{a+2uj,f} 

-B^pha(c7 + 2uj,){-4-(B-'fa^ + ph(5c7^+ 18cra;/+ 16^,^)}, (A17) 
ao = -2B^c7^uj,ia + 2uj,f - 'Xp^h^LOiia + uj,)ia + 20;/)"* 

- 2B^p^h^auj,(a + 2uj,f(5a + 6a;/) 

- AE^ph(j^Lo,{a + 2uji)(2a^ + lauji + 6ujj). (A 1 8) 

This equation should include the propagation modes. 

Eq. ( |A15l l is the biquadratic equation with respect to k, but includes unknown quantity w/. Figs. lAll show that the propagation 
mode becomes light waves in the short wavelength region, and damping rate -ujj is a monotonically increasing function of k. 
Note that what we want to know is whether the undamped shortest wavelength mode is light waves or Alfven waves, and we do 
not necessarily have to solve the biquadratic equation directly. 

For this reason, we substitute -2tt/ At into luj of Eq. ( IA15I ), and solve it with respect to k^: 

a 



=(/3i + V/52)//33 (A19) 
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Pi = B\B''fa^ -B^(T^{(T + 2a;/) 

- IB^pha^ila^ + 7crw/ + 6uj]) + p^h\(7 + lujif [-4(B*") + Iphia +2uj,f] 

-B^pha(a + 2Lu,)[-4(B"fa^ + ph(5<T^+ l^auj,+ 16Loj)] (A20) 
P2 = -8 p^h^uj,(a + lujif [B^a + 2ph(a + w/)] [B V + ph(a + 2uJi)f 
+ [-B (cr + 2w/) + 2p2/z2(cr + 2uJif{2{B''fa^ -ph(a+ 2uj,f } 
+ B'^a^iiB^fa^ - 2ph{2(j^+lauJi + 6uj])} 

+ B^pha{(T + 2uj,){A{B''fa^-phi5a^+\iauJi+ I6uj])}f, (A21) 
/33 = p^h^ia + 2L0i)[B^a + ph{a+2ui)], (A22) 

where Eq. dAlSI l has two solutions of k^, and we adopt the larger one since A; is a real number 

Substituting above into Eq. ( IA14b . one can obtain the desired characteristic velocity. Since what we need is appropriate 
MHD characteristic velocity, obtained velocity cannot be used as the characteristic velocity. However, if obtained velocity is not 
Alfven velocity, it shows that we should use light velocity for the characteristic velocity. 

This method requires some further explanation. 

First, note that the above method needs B^, B\ and At" in the comoving frame, and one should transform numerical data from 
Lab frame to comoving frame. 

Second, numerical experiments indicate that becomes for some range of k for some parameter, and of Eq. ( |A19l l 
becomes negative. In this case, we use speed of light as the characteristic velocity. 

Finally, Fig. lAll implies that -uji/a has some maximum value. This can be proved as follows. First, dividing Eq. jAlOl l by a'^, 
one obtains 

phui'* + i(B^ + 2ph)u? - [k^ph + (B^ + ph)]Lj^- i(B^ + ph)k^Lj + {B^fk^ = 0, (A24) 
where u = lu /a, and k = k/a. 

Fig. lAll implies that the maximum value of -w/ /cr is obtained in the limit of large k, and the propagation mode is light wave in 
this limit. For this reason, we substitute uj = k-idJ/ into Eq. (IA241 l. Then it reduces to 

-ip(-l +2uj/)P + [(B''f + B\-l +2uj/) + p(-l +5cj/-5w/2)]P 
+ i[-B^u},'i-2 + 3uj/) + 2puj/(l - 3uj/ + 2uj/^)]k 

-B\-1 +uj/)u;/^ + p(-l+uj/fLu,'^ = 0. (A25) 

Since this is in the limit of large k, what we have to consider is only the highest degree of k. Then, we set its coefficient equal to 
0, and it reduces to 

w/=^. (A26) 

This shows that -lu; / a becomes 1 /2 in the limit of large k/a, and we use the speed of light as the characteristic velocity when 
oj/ = -211 1 At is less than -g 12. 

MUSCL 

For the second-order scheme, one has to compute the cell boundary numerical flux using Riemann solver or metho d of char- 
acteri stics with left and right states obtained by using MUSCL. In this section, we explain MUSCL of Van Leer (Ivan Leer I 
[1979'). 

Since we need a second-order scheme, the left and right states of primitive variables Q are 

„«+l/2 _„«+l/2^2? 
^n+\ll _^«+l/2 5Q"i+\ 



^i""^ = f^"- [^(e;ViAL)-^(eti/2,«)], (B3) 



where follows from a predictor step 

At" 
2A^i 

where U is the conserved variables. In the above equation, 2"±i/2 '^^'^ be computed from Eqs. (IB 11 1 and (IB2l i by replacing Q"'^^/^ 
with Q". 

When one uses MUSCL, the SQt in Eqs. (IB lb and ( IB2I ) are computed as follows: 

min (2 1 AQm/2 \,\AQi\,2\ Afi,-! /2 1 ) sgn AQ/ 
(SQi)mono={ if sgnAa+i/2 = sgnAa = sgnAa_i/2, (B4) 
otherwise, 
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where 



Aa+i/2 = a+i-a, 



I/: 



(B5) 
(B6) 
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